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Abstract 

As a feasibility study for a scaling test we investigate the behavior of 
algorithms for dynamical fermions in the Nf = 2 Schrodinger functional at 
an intermediate volume of 1 fm . Simulations were performed using HMC 
with two pseudo-fermions and PHMC at lattice spacings of approximately 
0.1 and 0.07 fm. We show that some algorithmic problems are due to large 
cutoff-effects in the spectrum of the improved Wilson-Dirac operator and 
disappear at the smaller lattice spacing. The problems discussed here are 
not expected to be specific to the Schrodinger functional. 



1 Introduction 



1.1 Motivation 

For improved Wilson fermions it has long been established that in the quenched 
approximation cutoff-effects at a lattice spacing of a~0.1fm are tolerable and a 
continuum extrapolation can be started there. Recently more and more evidence 
has been accumulated that for dynamical improved Wilson fermions in a similar 
physical condition the cutoff-effects are much larger than expected. As an extreme 
example, for three flavors the existence of a phase transition in the /3-K-plane 
has been numerically conjectured and is interpreted as a lattice artifact [1]. A 
summary of large scaling violations in the two-flavor-theory is given in ref . [2] . 

In order to quantify those we are preparing a scaling test similar to what was 
done for the quenched case in ref. [3] . This will also serve as a benchmark for new 
actions. In the course of the scaling study we plan to calculate the axial current 
normalization constant Za(<7o) an d the axial current improvement constant ca(<7o) 
using the methods described in refs. [4] and [5], respectively. 

On our coarser lattices we encountered algorithmic difficulties in both the 
molecular dynamics integration of the Hybrid Monte Carlo (HMC) and the ef- 
ficient simulation of the canonical ensemble. We thus found it advantageous to 
deviate from importance sampling. Here we will discuss these problems and their 
link to cutoff-effects at the infrared end of the spectrum of the Dirac operator. 

This paper is organized as follows: In the remainder of Section^Jwe will briefly 
describe our setup and give a summary of the parameters of the simulations, which 
we will quote in the following sections. In Section |21 we describe in more detail 
the problems encountered at large coupling and also discuss methods to address 
these. In this context we study the behavior of the Polynomial Hybrid Monte 
Carlo (PHMC) algorithm [6,7] in this situation and find it a very useful tool for 
a detailed investigation of the properties of the small eigenvalues of the Dirac 
operator. 

Section |3] is devoted to a comparison of the spectrum with the quenched 
case at matched physical parameters and in Section 0] we present results from an 
exploratory run at a smaller lattice spacing. We close with a brief summary of 
our findings. 

1.2 Setup 

All our simulations were performed in the Schrodinger functional (SF) setup [8,9]. 
We use non-perturbatively improved Wilson fermions [10-13] and the plaquette 
gauge action. For Nf = 2 the clover coefficient csw at /3 = 5.2 has been set to the 
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value suggested in [12] and recently confirmed by the JLQCD Collaboration [13]. 
For the additional boundary-improvement coefficients needed in the SF we used 
the perturbative values for c t (2-loop) [14] and c t (1-loop) [15]. The axial current 
improvement constant ca is also set to its 1-loop value [16]. The first algorithm 
used is the HMC with two pseudo-fermion fields as proposed in ref. [17]. We 
want to note that the physical situation here is quite different from the one where 
this algorithm was previously applied and its performance tested by the ALPHA 
Collaboration [18]. In our planned scaling study we are interested in intermediate 
size physical volumes and lattice spacings between 0.1 and 0.05 fm. As mentioned 
above the second algorithm we employed is the PHMC, which we will discuss in 
some detail in Section |21 Apart from global sums all our calculations are carried 
out in single-precision arithmetics. 

1.3 Simulation parameters 

In Tabled we list the lattice sizes and bare parameters of our simulations. In all 
cases we have T = 9/4 L for the time extension T. The bare quark mass m is 
defined in the appendix of ref. [19]. In the algorithm column 'Ha' refers to HMC 
with two pseudo-fermion fields and 'P n ' stands for PHMC with a polynomial of 
degree n. The trajectory length is always equal to one and the molecular dynamics 
integration step-size is denoted by St. For each simulation we ran 16 independent 
replica to gain more statistics. Concerning the SF parameters we work with zero 
background field and periodic spatial boundary conditions (# = 0). 



run 


L/a 


P 


K 


csw 


Lm 


algo. 


A'traj St 


acc. 


I 


8 


5.2 


0.13550 


2.017 


0.205(10) 


H 2 


16-500 1/16 


91% 


II 


8 


5.2 


0.13515 


2.017 


0.307(9) 


H 2 


16-520 1/25 


97% 


III 


8 


5.2 


0.13515 


2.017 


0.314(8) 


Pl40 


16-500 1/26 


87% 


IV 


8 


5.2 


0.13550 


2.017 


0.195(7) 


Pl40 


16-400 1/25 


85% 


V 


8 


6.0 


0.13421 


1.769 


0.193(3) 




- quenched - 




VI 


12 


5.5 


0.13606 


1.751 


0.287(3) 


H 2 


16-240 1/20 


91% 


VII 


12 


6.26 


0.13495 


1.583 


0.295(3) 




- quenched - 





Table 1: Summary of simulation parameters. 
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2 Sampling problems on coarse lattices 

2.1 Instabilities in the molecular dynamics integration 

Algorithms making use of molecular dynamics (MD) require a numerical integra- 
tion of the equations of motion. Along a trajectory the Hamiltonian is then only 
conserved up to powers of the step-size 5t employed in the integration. Apart 
from these small deviations, under certain conditions the currently used integra- 
tion schemes can become unstable and produce very large Hamiltonian violations 
AH. For a more detailed discussion see ref. [20], where a connection between 
these instabilities and large driving forces in the MD is proposed in analogy to a 
harmonic oscillator model. In this model the integrator becomes unstable when 
the product of the force and the integration step-size exceeds a certain value. 

The reversibility of the numerical integration is needed to prove detailed bal- 
ance for these algorithms, which in turn implies that (e~ AH ) = 1. Here one should 
note that the average is taken over all proposed configurations (see ref. [18]). 
Therefore this quantity is also sensitive to those, which were rejected in the 
Metropolis step following the MD integration, i.e. trajectories resulting in a large 
value of AH. In a histogram of e~ AH these contribute to bins close to zero while 
the distribution is peaked around one. They can also lead to an unusual auto- 
correlation of this quantity, making the Monte Carlo error estimate difficult. 1 In 
particular this applies also to the integrated autocorrelation time of e~ AH itself. 
This is due to the long periods of rejection in the Metropolis step, which sometimes 
follow large AH values. 

Fig. ^ shows a histogram of e~ AH and also its integrated autocorrelation time 
from one of our simulations. In this data set there are several series of large AH 
values, during which the proposed configurations were rejected. In the distribution 
of e~ AH these lead to an additional peak close to zero. One also sees from the 
right-hand plot that e~ AH is noticeably autocorrelated only when a large number 
of proposals were rejected in a row. As argued above in these cases the error of 
Tint could be underestimated. These two effects might cause some concerns when 
using (e~ AH ) — 1 as an indicator for the absence of reversibility violations [18]. 

Spikes in AH have been observed by several collaborations using (improved) 
Wilson fermions in various setups (e.g. different gauge actions and volumes) at 
relatively large lattice spacings [12,20,22,23]. There these spikes have been traced 
back to large values of the driving force in the MD evolution and also their de- 
pendence on the quark mass has been investigated. 

1 All our data analysis is done using an explicit integration of the autocorrelation function as 
detailed in ref. [21]. This method also provides an estimate of the error of r m t- 
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Figure 1: e from run I. Left plot: histogrammed distribution. Right plot: 
our estimates of T mt in units of MD time separately for the 16 replica. In our 
normalization r int = 0.5 means no autocorrelation (dotted line). 



Here we want to clarify a point, which is essentially implied by the previous 
observations [22,24], namely the strong correlation between spikes in AH and 
small eigenvalues of the Dirac operator. 2 In this way we hope to be able to 
separate physical effects from cutoff-effects, i.e. the occurrence of unphysically 
small eigenvalues. In Fig.|2]we clearly see a long period of rejection (corresponding 
to the rightmost data point in Fig. caused by the presence of a very small 
eigenvalue. Although we did not measure them, this is expected to produce large 
fermionic contributions to the driving forces since they involve an inverse power 
of the Dirac operator. 

We found the observed average A m i n to be close to its tree-level estimate 
with Schrodinger functional boundary conditions [9]. However, the smallest X min 
is an order of magnitude below that and we therefore consider these eigenvalues 
unphysical and will later establish their nature as cutoff-effects. 

Finally, following the procedure of ref. [18], the absence of global reversibil- 
ity violations is explicitly verified even for trajectories resulting in large values of 
AH. Nevertheless our experience shows that the increased cost of using a smaller 
5t such that no long periods of rejection occur is more than compensated by the 
reduction in autocorrelation time of all observables. The reason is that already a 
small decrease of the integration step-size greatly reduces the Hamiltonian viola- 
tions. For example, repeating run I with a step-size of 1/20 instead of 1/16, the 
longest period of rejection was 4 (instead of 47) consecutive trajectories. 

2 Here and in the following we will always refer to the eigenvalues of the square of the Hermi- 
tian even-odd preconditioned Dirac operator Q 2 in the Schrodinger functional. For its precise 
definition see rcf. [18]. 
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Figure 2: Monte Carlo history for one replicum of run I with a long period of 
rejection. Configurations where the new proposal was rejected are marked by 
a dot. At £md= 340 the algorithm gets stuck with a configuration carrying an 
exceptionally low smallest eigenvalue A min of Q 2 . 



2.2 MC estimates of fermionic observables 

We concluded in the previous section that unphysically small eigenvalues of Q 2 
produce algorithmic problems only on a practical and not on a theoretical level. 
But apart from slowing down the algorithm these small eigenvalues also cause 
problems in the MC evaluation of fermionic Green's functions. 

Consider the Schrodinger functional correlation function j\ as defined in 
ref. [19]. It is the correlation between pseudo-scalar composite fields at the first 
and last time-slice, respectively. We will denote its value on a given gauge field 
configuration by [0i]f- Fig- El shows the MC history of the normalized [4>\]f for 
the 16 replica of run II. Here ^md refers to the molecular dynamics time for each 
replicum. While on this scale the bulk of the data are below one and hence not 
visible there are several peaks, which have a big influence on the mean value. 
These spikes also affect the error estimate cr(fi) through both the variance and 
the integrated autocorrelation time [21]. For statistically accessible quantities the 
error should approach a 1 / v^md behavior in the limit £md — >• oo. In this respect 
we found f\ and all other fermionic correlation functions we considered to be very 
hard to measure. Even when using 16 replica, this asymptotic behavior does not 
set in after t MD ~500. 

The reason is the rare occurrence of very large values of [4>i]f, which appear 
to be correlated with small eigenvalues of Q 2 . However, this effect is washed out 



6 




100 200 300 400 500 



Figure 3: Monte Carlo history for the iV r ep = 16 replica of run II showing the 
normalized [0i]f- 

by using several replica. We therefore show in Fig. 0] the MC history of [4>i]f, 
A m j n and our error estimate for fi for one replicum of run II with such a spike 
in [0i]f- Indeed, for each spike in [0i]f the smallest eigenvalue drops below its 
average. That the converse is not true could be ascribed to a lack of overlap 
of the eigenvector corresponding to A min with the source needed to compute the 
quark propagator. Quantitatively, for the correlation between [0i] F and A min we 
measure a value of C*[</ )1 ] F ,A min = —0.33(4) if we use all replica and —0.46(6) from 
the replicum shown in Fig. 0] alone. Here we used as a definition of the correlation 
Ca,b between two observables A and B 



C A , B = , ( AB )~( A )( B ) sothat -1<C AB <1. (2.1) 



Even though in the limit of infinite statistics configurations carrying very 
small eigenvalues are given the correct weight, depending on the algorithm this 
might be badly approximated for a typical ensemble size. Similar arguments 
referring in particular to the HMC algorithm motivated the introduction of the 
Polynomial Hybrid Monte Carlo (PHMC) algorithm in refs. [7]. 

Hence the difficulty in measuring fermionic correlation functions might be an 
efficiency problem related to the choice of the algorithm. To check this conjecture 
we employ a second algorithm and compare ensembles generated by HMC (with 
two pseudo-fermion fields) with PHMC ensembles. Indeed, PHMC can be tuned 
in such a way that it enhances the occurrence of configurations carrying small 
eigenvalues, thus resulting in a better sampling of this region of configuration 
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Figure 4: Normalized [</>i]f and smallest eigenvalue from one "sick" replicum of 
run II. Evidently the spike in [0i]f is dominating the statistical error cr(fi). 



space. A reweighting step is introduced to render the algorithm exact. As a 
preparation for the following discussions we want to recall some properties and 
introduce the notations concerning the PHMC 



2. 2. 1 The PHMC algorithm 

One of the main ideas of the PHMC algorithm is to deliberately move away from 
importance sampling by using an approximation to the fermionic part of the lattice 
QCD action. More precisely, in an HMC algorithm the inverse of Q 2 is replaced 
by a polynomial P Tlt<i (Q 2 ) of degree n. Here P n>€ {x) approximates 1/x in the range 
e < x < 1. As a consequence this algorithm stochastically implements the weight 
dU det P~l(Q 2 )e~ Scj , whereas standard HMC generates ensembles according to 
dU det Q 2 e~ Sa with S g being the gauge part of the action and U the gauge link 
configuration. Denoting averages over the PHMC ensemble by (. . .)p, the correct 
sample average of an observable (O) can then be written as 

W 

(O) = (Olj) p , where u = — — , (2.2) 

\W)p 

and we introduce the reweighting factor W as a (partially) 3 stochastic estimate of 
det{Q 2 P ntt (Q 2 )} . When using Chebyshev polynomials the relative approximation 
error for e < x < 1 is bounded by 5 ~ 2 exp(— 2y/en). 

3 Through the separate treatment of the lowest eigenvalues of Q 2 the infrared part of W is 
evaluated exactly. 
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Figure 5: Three different Chebyshev polynomials approximating 1/x, all with 
5 = 0.001. The right-hand plot shows the relative deviation from 1/x as a function 
of x. There the points (e, 5) are marked by dots. 



To give an impression of the role of e and 5 we plot in Fig.EJa set of polynomi- 
als P n ^ e {x) for typical (in our simulations) values of these parameters and compare 
them with 1/x in the region of small x. Depending on the smallest eigenvalue of 
Q 2 the parameters e and n have to be tuned such that the reweighting factor does 
not fluctuate too much. The authors of ref. [7] suggested to take e of the same 
order as (A m i n ) and in practice used e ~ 2(A m j n ) and 5 < 0.01. 

Recalling that PHMC replaces detQ 2 in the HMC weight with det P~l(Q 2 ) 
and observing from Fig.EJthat P n>e (x) is smaller than 1/x for x < e, the aforemen- 
tioned property of enhancing the occurrence of small eigenvalues is evident. At 
this point we would like to note that the fermionic contribution to the driving force 
in the PHMC is bounded from above since P n>e (x) is finite even at x — 0. In this 
way the polynomial provides a regularized inversion of Q 2 , thus also addressing 
the problems mentioned in Section 12.11 



2.2.2 HMC vs. PHMC 

Coming back to the comparison of samples from HMC and PHMC, we repeated 
run II with PHMC using a polynomial of degree 140 and e = 6 -10~ 4 , resulting in 
5 ~ 0.002. The ratio e/ (A m i n ) turned out to be around 2.7. In Fig. |U1 we plot for 
this run the MC history of [0i]p and of [0i]f ■ which enters into eq. (|2.2j) if we 
consider O = [0i]f, i-e. 

h = <Nf> = (Nf • w)p = {[ ^}1' W)P ■ (2.3) 
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Figure 6: Monte Carlo history for the 16 replica of run III showing the correlation 
function [0i]f and the product [0i]f ■ where uj is the normalized reweighting 
factor. Our error estimate of j\ shows the expected scaling behavior as soon as 
the run is long enough for a reliable extraction of r^f 

We first observe that apart from removing the largest spikes the inclusion of the 
reweighting factor does not seem to significantly change the relative fluctuations. 
This means that the parameters of the polynomial have been chosen properly. 
Events where [0i]f assumes a value O(10) times larger than fi are no longer 
isolated as in Fig.|3]but happen frequently, which means that the PHMC algorithm 
can more easily explore the associated regions in configuration space. This is what 
allows a reliable error estimate as shown in the upper part of Fig. El i.e. with 16 
replica the asymptotic behavior of the error sets in after £md — 100. 

The advantage of using PHMC instead of HMC can be clearly seen by con- 
sidering the spread of ct(/i)v / ^md among the replica. We analyzed this quantity 
in extensions of runs II and III. The result is shown in figure Fig. where the 
shaded areas represent the range of values covered by the 16 replica as a function 
of the MD time. In the limit of infinite statistics all replica should converge to 
the same value, which need not be the same for the two algorithms because of 
reweighting and different autocorrelation times. We see that the spread for the 
HMC data is more than twice as large as for PHMC, i.e. the error on j\ is sig- 
nificantly harder to estimate with HMC. What we are suggesting here is that the 
algorithm should be chosen depending on the type of observables and the param- 
eter values. From our experience we conclude that PHMC sampling might just be 
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Figure 7: Monte Carlo history of o"(/i)v^md from extensions of runs II and III. 
For the two algorithms we show the ranges covered by the 16 replica. 

more effective than HMC when computing fermionic quantities that are sensitive 
to small eigenvalues. 

To gain some more insight into the difference in sampling we consider the 
distribution of A m j n since this is where we expect the largest effect. The distri- 
butions are analyzed by treating Abm = Xbin(A m in) as a primary observable. Here 
Xbm denotes the characteristic function of each given bin in the histogram. We 
then perform our normal error analysis for (A^n), where eq. (|2.2j) has to be used 
if it is a PHMC sample. For comparison (Abm)p is also analyzed in this case. 

The histograms in the upper part of Fig. |H] compare the results from 200 
independent measurements produced by HMC and PHMC (runs II and III, re- 
spectively). As expected the distributions agree within errors. For the PHMC 
run we also plot the unreweighted histogram, i.e. (Abin)p- Here we again con- 
firm that with the parameters we chose for the polynomial the PHMC produces 
more configurations with small eigenvalues than HMC. As a consequence of the 
reweighting the errors at the infrared end of the spectrum should be smaller for 
the PHMC data. This is explicitly verified in the lower part of the plot where we 
show the ratio of the errors on (Abi n ) from the two algorithms. The three symbols 
refer to different bin sizes. The advantage in using PHMC to sample this part of 
the spectrum is significant and we will make use of this in the following discussion. 
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Figure 8: Upper panel: histograms of A m i n , i.e. (Abm) vs. 'bin', from runs II and 
III. For run III we also show (Abm)p- Lower panel (from the same data): ratio of 
the error on (A^n) from HMC to that from PHMC 



3 Comparison to the quenched case 

In the previous section we studied various problems related to the occurrence of 
small eigenvalues. All the data presented there were produced at bare parame- 
ter values, which correspond to relatively large quark masses and small volumes. 
These small eigenvalues might therefore have a different nature from the "physi- 
cal" ones expected to show up in large volumes and/or close to the chiral limit. 
Here and in the next section we will establish them as cutoff-effects. 

To this end we made an additional simulation at the parameters of run II 
and calculated the ten lowest-lying eigenvalues Aj, i = 1 ... 10. In Fig. |H1 the 
smallest eigenvalue, Ai, is denoted by an open symbol. It seems that while A2 
through A10 form a rather compact band, the lowest eigenvalue fluctuates to very 
small values quite independently of the others. It is expected and has been shown 
numerically [25] that the spectrum of the Dirac operator depends quite strongly 
on the bare gauge coupling. A well-defined lower bound should be recovered close 
to the continuum limit only. Therefore we take the strong fluctuations of A min 
as an indication for the presence of large cutoff-effects. Here we should point 
out that the eigenvalues of the Dirac operator are not on-shell quantities. Hence 
the Symanzik improvement programme does not necessarily reduce cutoff-effects 
here. Quenched experience even suggests that the opposite might be true [26]. 

The occurrence of small eigenvalues at these bare parameters poses a some- 
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Figure 9: Monte Carlo history of the ten lowest eigenvalues at the parameters of 
run II. The open symbols denote A m i n . 

what unexpected problem in dynamical simulations. Comparing the quenched 
situation to the Nf = 2 dynamical case, the naive expectation is that at fixed 
bare parameters the probability of finding configurations with small eigenvalues 
should be reduced by the determinant. To us the more relevant question seems 
to be whether small eigenvalues are suppressed in a situation where the physical 
parameters (e.g. volume and pseudo-scalar mass) are kept constant. 

Using the quenched data from ref. [27] and the dynamical data from refs. [28] 
and [29] (where an estimate of r^/a — 5.21(6) for [3 = 5.2 can be found) we chose 
the parameters of the quenched run V such that the lattice spacing and the (large 
volume) pseudo-scalar mass are matched to run IV. This was found to occur at 
almost equal bare current quark mass (see Lm in Table [1}. In Fig. ^] we compare 
the distributions of A m j n for these two runs. Two comments are in order here: 

• For the dynamical run the mean value is shifted up from 1.44(1) ■ 10 -4 to 
1.72(5) • 10~ 4 . This agrees with the naive expectation but in a physically 
matched comparison it is a non-trivial observation. 

• The distribution itself is significantly broader compared to the quenched 
case and in particular it is falling off more slowly towards zero. This means 
that even though (A min ) is larger for Nf = 2 the probability of finding very 
small eigenvalues is enhanced. 

The second point, i.e. that the lower bound of A min is less well-defined, seems 
to imply that at a lattice spacing of a ~ 0.1 fm the cutoff-effects are much larger 
in the Nf = 2 case. To substantiate this we will compare the distribution of A min 
to that from a run at finer lattice spacing and matched physical parameters. 
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Figure 10: The histograms of A m j n from run IV (Nf = 2) and run V (Nf = 0). 
Despite its higher mean value the dynamical data show a much larger probability 
of finding very small eigenvalues. 



4 Finer lattices 



Apart from cutoff-effects, in the massless theory the Schrodinger functional cou- 
pling g 2 is a function of the box size L only [8,9]. We measured it on a small 
lattice of extension L/a = 4 at = 5.2, obtaining a value of g 2 = 3.7(1). We then 
extrapolated to this value the L/a = 6 data used in ref. [30] as a function of 0. Our 
result from the matching is that for the two-flavor theory a bare gauge coupling 
of /? = 5.5 roughly corresponds to a lattice spacing, which is 1.5 times smaller than 
at = 5.2. 

Hoping that the algorithmic difficulties arising from cutoff-effects would be 
much smaller in this situation, we simulated a 12 3 x27 lattice at this value of 
(run VI) using the HMC algorithm. With the k we chose (and ignoring the change 
in renormalization factors) the bare quark mass Lm is roughly matched to the 
heavier runs at = 5.2. We therefore compare run VI with run III. 

Normally, a constant acceptance requires a decrease of the MD integration 
step-size if ones goes to finer lattices at fixed physical conditions. This argument 
is based on the scaling of the small eigenvalues, which influence the MD driving 
force. We found that (A m i n ) in run VI is a factor two smaller than in run III. 
Nevertheless, at = 5.5 the step-size necessary for a certain (~ 90%) acceptance 
is roughly the same as at = 5.2. This indicates that the value of 5t we had to 
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Figure 11: Normalized distributions of A m i n from runs III (/3 = 5.2) and VI (j3 = 
5.5). While the data from the coarse lattice stretch almost to zero, the (3 = 5.5 
data seem to have a more well-defined lower bound. 

use in the HMC runs at (3 — 5.2 was dictated by the occurrence of extremely small 
eigenvalues rather than by the average smallest eigenvalue. In addition, where in 
run I at the same average acceptance a maximum of 47 proposals were rejected 
in a row, the maximum for run VI is 4 trajectories. For this reason e~ AH shows 
no autocorrelation and its distribution is well separated from zero. 

Concerning fermionic observables, we have not observed spikes and hence 
expect the error to scale properly. However, for an accurate estimate of the error 
on e.g. fi our present statistics is not yet sufficient. 4 

The reason for these effects is the change in the distribution of A m ; n . To com- 
pensate for the different lattice spacing, Fig. ITT1 compares A min /(A min ) from runs 
III and VI. One can clearly see that at the finer lattice spacing the probability of 
finding a smallest eigenvalue less than half its average is greatly reduced compared 
to /3 = 5.2. The width of the distribution is smaller in this case and in particular 
the spectrum is now clearly separated from zero. Quantitatively, the normalized 
variance of A min is reduced from 0.18(1) to 0.13(2). 

This comparison explicitly shows that the long tail of the eigenvalue distribu- 
tion we observed at a ~ 0.1 fm, and which caused the problems we have discussed, 
is a cutoff-effect. Matching also run VI to a quenched simulation (run VII), we 
again found an upward shift of (A m i n ) for the dynamical case. In addition, at this 

4 Ratios of correlators relevant for physical applications are easier to estimate. 
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Figure 12: The histograms of A m i n from run VI (Nf — 2) and run VII (Nf = 0). At 
this finer lattice spacing the lower end of the spectrum appears to be similar in 
the quenched and the dynamical case. 



finer lattice spacing, the tails of the distributions of A m i n look already very similar 
to each other as shown in Fig. 



5 Conclusions 

At a lattice cutoff of approximately 2 GeV we have studied the behavior and per- 
formance of HMC-type algorithms in an intermediate size volume of lfm 4 . We 
discussed problems related to the occurrence of small eigenvalues in two-flavor 
dynamical simulations with improved Wilson fermions. We found these small 
eigenvalues to be responsible for large Hamiltonian violations in the molecular 
dynamics. Even for integration step-sizes such that the acceptance is 80 ~ 90% 
those can still cause long periods of rejection, thus degrading algorithmic perfor- 
mance. However, in spite of employing only single-precision arithmetics we never 
observed reversibility violations. 

In addition, those eigenvalues make the estimate of fermionic quantities very 
difficult. The naive intuition is that the fermionic determinant should suppress 
small eigenvalues compared to the quenched case. Through a direct comparison 
at matched physical parameters we indeed verified that (A min ) is larger with two 
dynamical flavors. On the other hand there is no obvious expectation for the tail of 
the distribution and we observed that it extends further towards zero than in the 
quenched case. Given the infrared cutoff induced by the Schrodinger functional 
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boundary conditions and the quark mass we interpret this as a lattice artifact. 
We were able to confirm this picture with a simulation at finer lattice spacing, 
where the spectrum turned out to have a much sharper lower bound. 

In our study we found that the PHMC algorithm is more efficient than HMC 
(with two pseudo-fermions) in incorporating the contribution to the path integral 
of configurations carrying small eigenvalues. In other words, the distortion of 
the spectrum by cutoff-effects actually makes it advantageous to deviate from 
importance sampling. Also without such special problems we found PHMC at 
least comparable in performance to HMC (in our implementations). 

We want to emphasize that the problems discussed here do not occur only 
in the Schrodinger functional setup. Without this infrared regulator they are ex- 
pected to show up already at larger quark masses. 
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